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Abstract 

This paper proposes a methodology to calculate both the first and second derivatives of a vector function 
of one variable in a single computation step. The method is based on the application of the dual numbers. 
It has been implemented in Fortran language, a module which contains the dual version of elementary 
functions as well as more complex functions, which are common in the field of rotational kinematics. Our 
basic numerical entity has three elements, then, for a given vector function / : M — > M™, its dual version will 
have the form / : — >■ M^™. As a study case, the proposed methodology is used to calculate the velocity 
and acceleration of a point moving on the coupler-point curve generated by a spherical four-bar mechanism. 
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1. Introduction 

The calculation of velocity and acceleration is highly needed in the fields of physics and engineering. 
In some cases, there is no possibility to obtain them analytically. In other cases, the analytical result is 
quite complicated to deal with; in both situations a numerical treatment is desirable. Regarding the field 
of mechanisms, the calculation of the first and second derivative, lets the designer to include the velocity 
and acceleration in the synthesis of balanced mechanisms, dwell mechanisms, etc., [1, 2, 3, 4, 5, 6, 7]. 

Usually, the process of calculating a derivative is not difficult. However, for the case of a spherical 
mechanism, obtaining first and second derivatives of the position vector for the coupler point is not simple. 
Even when such derivatives can be explicitly obtained, the resulting expressions could be of great complexity 
and useless for practical purposes. So, an alternative solution is to numerically calculate such derivatives. 
Nevertheless, traditional methods for calculating numerical derivatives (finite-difi'erences) are subject to 
both truncation and subtractive cancellation errors. There is no doubt that methods without such problems, 
would be of great importance in order to obtain more precise numerical derivatives. Such methods have 
been developed using dual numbers. 

Analogous to the definition of a complex number z =^ a + ib with a and b being real numbers and 
i"^ — —1, a dual number is defined as f = a + eb where, a and b are real numbers but = 0. Although the 
first applications of the dual number theory to the development of mechanical engineering dates back to 
late XIX century [8, 9], and that the algebra of dual numbers was also developed in the late XIX century 
[10], applications of dual numbers to numerical calculation of derivatives are relatively recent [11, 12, 13]. 
Since then, a big amount of work regarding applications of dual numbers has been made. They are used 
to describe finite displacements of rigid and deformable bodies, for the analytical treatment in kinematic 
and dynamics of spatial mechanisms, for the study of the kinematics, dynamics and calibration of open- 
chain robot manipulators, for the study of computers graphics, etc. In references [13, 14, 15] there is some 
bibliography about such studies. 
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Due to the great applicability of the dual numbers, it is important to develop algorithms implementing 
their algebra, this has been done in [14]. It is worthwhile to mention reference [15], where some linear 
dual algebra algorithms are presented. Regarding applications of dual numbers to compute the numerical 
derivative of functions, we can cite [11, 12, 13]. In [11], dual numbers are used to calculate first order 
derivatives and in [12] and [13], they are used to calculate second order derivatives by using the operator 
overloading method. 

The aim of this paper is twofold. First, it develops a methodology based on the dual numbers, where 
first and second order derivatives are obtained in a single computation step. Second, it obtains the velocity 
and acceleration of a point moving on the coupler-point curve generated by a spherical four-bar mechanism. 

Once the methodology to develop elementary functions in its dual version including its second derivative 
has been presented, it is straightforward to develop more sophisticated functions, such as rotations and 
numerical solutions to equations. Moreover, as the use of functions in programming languages as C, CH — h 
and Fortran is quite intuitive and easy to codify, we have created a Fortran module where some functions 
in its dual version, as well as, common functions to the field of rotational kinematics, are provided. Such a 
module along with an example of its usage is presented in the appendix. 

The rest of the paper is organized as follows. Section 2, presents the process of calculating the numerical 
derivative using dual numbers. Section 3, presents the implementation in Fortran language of the dual 
versions of scalar elementary functions, as well as, more complicated expressions, such as vector and matrix 
functions. In section 4, we show how the velocity and acceleration of a point moving on the coupler- 
point curve generated by a spherical four-bar mechanism are obtained. Finally, section 5, presents the 
conclusions. 

2. Dual numbers and derivatives 

A dual number f is a number of the form 

f^a + eb, (1) 

where a (the real part) and b (the dual part) are real numbers and = 0. As in the case of complex 
numbers, there is an isomorphism^ between dual numbers and the real vector space M^. So a dual number 
f can be defined as ordered pairs 

r^{x,y}. (2) 

The algebraic rules for dual numbers can be found elsewhere in the literature, see for example [10, 16, 15]. 
Below we present how can the dual numbers be used to calculate the derivatives of a function. 

2.1. First order derivative 

Let us consider the Taylor series expansion (3) of a function / : M — > M, around the point x, where h.o.t. 
stands for higher order terms. 

fix + h) = fix) + f'ix)h + ^^h^ + ■■■ + h.o.t. (3) 

Now, let us consider the dual number x = x + e, i.e., a dual number where the coefficient of the nilpotent 
e is equal to one. Substituting x in (3), we obtain (4). Then, if instead of computing over the reals, we 
compute over the dual numbers, we came up with a dual function / = fix + e) whose real term is the 
original function fix) and the coefficient of e is its derivative fix). Thus 

/(.T) = /(.T + e) = /(.T)+/'(x)e, (4) 



^Undcr the ordinary addition and scalar multiplication — multiplication by a real number, the transformation T(a + be) = 
{a, b} is linear and bijective. 
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or in the same notation of (2): 

= (5) 

Let us exemplify the procedure for computing the dual function, by using the sinusoidal function. Let 
f{x) — sin(a;) be the function for which the dual function will be obtained and let g{x) = {gi{x), g2{x)} be 
a dual function with gi{x) — g{x) and 52(2;) = g'i^). Then, from (5) we obtain (6). 

sm{g{x)) = {sm{gi{x)) ,cos{gi{x)) g2ix)} . (6) 

2.2. Second order derivative 

The second order derivative can be obtained by applying (5) to f'{x), that is to say, by writing the dual 
version of the function f'{x): 

f'{x)^{f'ix)J"ix)}. (7) 

The relevant information can be allocated in a vector of three components. Such a vector will have 
the information of the function, its first and its second derivative, so we could speak of an extended dual 
function (to differentiate it from the common dual function which only have two components). We will use 
the notation 

fix)^{f{x),.nx)J"{x)}, (8) 

to represent the extended dual version of the original function /. The identity function in its extended dual 
version i, has three components and it is given by (9). 

x^{x,l,0}. (9) 

Let us exemplify the proposed methodology, using the sinusoidal function. Let g{x) = {gi(x), g2{x), g^lx)} 
be an extended dual function, where gi{x) = g{x), g2{x) = g'{x), gsix) = g"{x). The sinusoidal function 
in its extended dual version is given by (10), where the arguments of the functions are not written to avoid 
unnecessary notation. 

sm(^) ={sin(gi), cos{gi) g2, -sv[i{gi)gl + cos(.gi) 53}. (10) 

Notice that no matter how complicated the function g, Eq. (10) ensures that the chain rule will be 
successfully applied. Thus by writing all the functions in its extended dual version the derivatives will be 
obtained without the need of the traditional methods of finite-differences. 



3. Fortran implementation of the extended dual functions 

The extended dual functions are defined as arrays of 3m components and since we are interested in 
functions / : M — ?> M™, a function g : K'^ — >■ K^™ is required to write the extended dual version of the 
function /. For example, in the case of the sinusoidal function and considering Fortran programming 
language, the code results as follows: 

module sphmodual 
contains 

function sindual(x) result (f_result) 
implicit none 

double precision :: f_result(3) 
double precision, intent(in) :: x(3) 
f_result=[sin(x(l)) ,cos(x(l))*x(2) ,& 
-sin(x(l))*x(2)**2 + cos(x(l))*x(3)] 
return 

end function sindual 
end module sphmodual 
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For coding purposes we will use f dual instead of /, so in the above code sindual means sin. A simple 
program calculating the function f{x) = sin(sin(a;)) for x = 1.1 in its extended dual version is: 

progrcun sind 
use sphmodual 
implicit none 

double precision :: angd(3) , res(3) 

angd=[l.ldO,ldO,OdO] 

res = sindual (sindual (angd)) 

print*, res(l) ,res(2) ,res(3) 

end progrcun sind 

After compiling and executing the above program, the result is: 
0.777831 0.285073 -0.720138 



The first component of the output vector is /, the second one is /' and the third one is /", all evaluated 
at the argument x ~ 1.1. The implementation of all the other elementary functions is straightforward from 
this example. 

3.1. Extended dual version of more complicated objects 

It has been shown how elementary functions can be dualized, but there are more complicated objects such 
as, for example, the vector product, whose extended dual version is required for applications of rotational 
kinematics. One approach to handle vector functions is to work by components. Let us exemplify for the 
cross product. The j-th component of the cross product of two vectors a and b is given in (11), where a„ 
and bn are the n-th component of vectors a and b respectively and Eijk is the Levi-Civitta tensor. Notice 
that summation over repeated indices is assumed. 

(a X b)j = £jjfeaj6fc. (11) 

As the derivative operator is a linear operator, only the extended dual version of the multiplication has 
to be obtained, but the addition has not. Below we show the Fortran code for the extended dual version of 
the cross product. 

function crossdual(xd,yd) result (f_result) 
implicit none 

double precision :: f _result(3,3) 

double precision, intent(in) :: xd(3,3) ,yd(3,3) 

integer : : j ,k 



f_result = OdO 



do j=l,3 

do k=l,3 
f _result(l , 
f _result(2, 
f _result (3 , 

end do 
end do 



f _result (1 , 
f _result(2, 
f _result(3. 



) + levied ,j ,k)*proddual(xd(j , 
) + levic(2, j ,k)*proddual(xd(j , 
) + levic(3, j ,k)*proddual(xd(j , 



),yd(k, 
),yd(k, 
),yd(k, 



)) 
)) 
)) 



return 

end function crossdual 
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In the levic function, the Levi-Civita tensor is coded. In the proddual function, the extended dual version 
of the scalar product is coded. The code for both functions is included in the appendix. 

The result is a 3 x 3 matrix where the first column is the real cross product, the second column is the 
first order derivative of the cross product and the third column is the second order derivative of the cross 
product. 

For example, let v(6') = {cos(6'), sin(0), 6*^} and w(6') = |c , 6'cos(0), sin(6')| be two vectors. A 
program that calculates its extended dual cross product at 6* = 1.1 is as follows: 

progrsmi crossp 
use sphmodual 
implicit none 

double precision :: aiig(3) ,xcl(3,3) ,& 
xc2(3,3) ,thr(3) ,crsp(3,3) 

thr = [3d0,0d0,0d0] 
ang = [I.ld0,ld0,0d0] 



xcld,:) = cosdual(cing) 
xcl(2,:) = sindual(cLng) 
xcl(3,:) = powdual(cLng,thr) 

xc2(l,:) = expdual(-proddual(cLng,aiig)) 
xc2(2,:) = proddual (cosdual (ang) , ang) 
xc2(3,:) = sindual(cing) 

crsp = crossdual(xcl,xc2) 
print* , crsp( : , 1) 
print* , crsp ( : , 2) 
print* , crsp ( : , 3) 
end progremi crossp 

Extended dual version of other mathematical objects like dot product, norms, matrix multiplications, 
rotation matrices, etc., can be implemented in a similar way. 



4. First and second derivative in the spherical 4R mechanism 

This section presents the application of the proposed methodology to compute the first and second order 
derivatives of some useful functions in the synthesis of mechanisms. In particular, the examples concern 
the spherical four bar mechanism, which is shown in Fig. 1. A detailed procedure in order to obtain the 
coupler-point curve can be found in [17]. Here we reproduce the essential formulas in order to make the 
paper self-contained. Notice that the notation has been slightly changed in order to have a clear code for 
the programming functions. 

4-.1. Derivative of the output angle 

Let us consider the spherical four bar mechanism shown in Fig. 1. Without loss of generality we assume 
a unit sphere. Such a mechanism is in its assembly configuration, hence the input angle is equal to Oq. The 
vectors x^; k — 1, 2, 3, 4 represent the initial position vectors for the joint k. We will consider the input link 
as the geodesic connecting the points xi and X2, the coupler link as the geodesic connecting the points X2 
and X3, the output link as the geodesic connecting the points X3 and X4, finally the fixed link will be the 
geodesic connecting X4 and xi. Since we are considering a unit sphere, ai, 012, as, will be the lengths of 
the input link, the coupler link, the output link and the fixed link respectively . Once the input link starts 
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its rotation, for example of an angle 0, its new position will be + 0. Similarly the position of the output 
link changes by the amount (j){9;'}i) where the symbol x represents the dependence on (xi, X2, X3, X4). In 
what follows the dependence of </> on x will not be written to avoid unnecessary notation, however it will be 
written in those vectors with such a dependence. 

Let r2{0; xi, X2) and r3(0(0); x) denote the final position for the extreme of the input hnk and the output 
link respectively. We have 

r2(0;xi,X2) =R(f?,Xi)x2 (12) 
r3(0(f?);x) =R(0(0),X4)X3, (13) 

where R('i?, x) is the active rotation matrix of angle 1} in direction of the unit vector x. 

Since the coupler link must have a constant length, the angle (f>{0) can be obtained from 

r2(6';xi,X2) • r3((/)(6'); x) = X2 • X3 = constant. (14) 
A closed-form solution to Eq. (14) can be found in [18, 19], in our notation, it is given by: 

0=$o-2tan-M 1, (15) 

where 

A = sinai sin 03 sm{6 + Qq) 

B = cos Qfi sin a3 sin — sin ai sin 0^3 cos cos{9 + Qo) 
C = sinai cos 0:3 sinQ;4 cos(d + 0o) + 
cos ai cos a3 cos — cos a2 ■ 

Although we use Eq. (15) to obtain the coupler-point curve generate by the mechanism, in this section, 
we are interested in showing a numerical approach, where the derivative of (j) with respect to can be 
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Figure 2: Vectors involved to obtain Tgon- The dependence on 9 and x is not written 



found without the analytical knowledge of the angle. This can be obtained by computing (16), where 

In order to obtain a numerical value for Eq. (16) the angle is obtained by numerically solving Eq. 
(14). Clearly the use of dual numbers is an advantage here, since by implementing the extended dual 
version of Eq. (14) and applying the numerical method of solution also in the context of the extended dual 
functions, the first and second derivatives of Eq. (14) with respect to 9 are automatically obtained, along 
with the solution for 4>. Notice that even when Eq. (16) is a formal solution to the problem at hand, yet 
the derivatives still need to be calculated. 

4-. 2. Velocity and acceleration of the coupler point 

We are interested in calculating the velocity and the acceleration of a point moving on the coupler-point 
curve. Given a mechanism, the only independent (excepting the time dependence) variable is the rotation 
angle of the input link, 9. Denoting as Tgon the position vector of the coupler point (see Fig. 2), we have 



(17) 



" _ h2 ^^(fgcn) -A C^(rgcn) „s 



for the velocity and the acceleration of the rgon vector respectively. Since one can, in principle, control the 
angular velocity of the input link, the problem is reduced to calculate d{rgcn)/d9 and d^{rgcn)/d9^. 

In order to obtain the parametric equation for the coupler-point curve, it is necessary to find the position 
vector rcp{9 ,v]ii) of an arbitrary point on the coupler link (see Fig. 2). This can be done by rotating the 
vector r2(^;xi,X2) by an angle u in direction of the unit vector 1123, which is orthogonal to r2(^?;xi,X2) 
and r3(0(6');x). Thus 

Vcp{9, v; x) = Ii{v, n23)r2(6'; xi, X2), (19) 

where 

^ ^ r2(6';xi,X2) X r3(0(6>);x) 

||r2(0;xi,X2)xr3(<^(0);5?)ir ^ ' 
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Table 1: Parameters of the Spherical 4R mechanism; the parameters (pi and rji are the azimuthal and polar 
angles for the point 
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92 




0.48867 


0.23066 


0.47437 


0.000009 


0.38828 


0.19646 


'■Pi 












0.97780 


1.57081 


1.46619 


0.66128 


1.34474 





Table 2: Components of the velocity vector, Ygcn{Q) (we have taken = 1 and the units are in the SI system) 



e 








0.00000 


-0.14254 


-0.06886 


0.59467 


0.62832 


-0.28007 


-0.18549 


0.35545 


1.25664 


-0.17826 


-0.23972 


0.05445 


1.88496 


-0.02698 


-0.20190 


-0.13156 


2.51327 


0.10680 


-0.11271 


-0.26667 


3.14159 


0.15217 


-0.00107 


-0.36591 


3.76991 


0.12510 


0.13205 


-0.36656 


4.39823 


0.09935 


0.25463 


-0.23445 


5.02655 


0.08944 


0.28060 


0.01395 


5.65487 


0.05510 


0.14224 


0.34717 



Using Eq. (19), the position vector of the coupler point Tgon is obtained by rotating the vector Vcp{9, /?+7; x) 
an angle 7r/2 in direction of the vector rcp{0, f3;x.): 

rgen(e, /?, 7; x) = R [71/2, rep(0, /?; x)) r,p{0, /3 + 7; x). (21) 

So, writing the Eq. (21) in its extended dual version, we could automatically get the parametric equation 
for the coupler-point curve, as well as, its first and second derivatives with respect to 9. It is worthwhile to 
mention that all the necessary functions to obtain rgcn in its extended dual form are coded in the module 
(see Appendix B for an example of its usage). 

As a practical example, let us consider the spherical four bar mechanism whose parameters are shown 
in Table 1, with the di parameter corresponding to the starting point (0.85737,-0.18481,0.48037). This 
mechanism was synthesized for the path generation task in [17]. The desired points were first presented in 
[20] and normalized in [21]. 

Tables 2 and 3, show numerical results for velocity and acceleration, respectively, of a point moving on 
the coupler-point curve as function of the input angle 9. 

5. Conclusions 

Velocity and acceleration of a point moving on the coupler-point curve generated by a spherical 4R 
mechanism, are obtained by using dual numbers. Although it could be possible to obtain an analytical 
expression for the position vector of the coupler point, the complexity of the expression does not allow an 
efficient method to obtain their derivatives, but, if the vector is written in its extended dual version, as it 
is proposed in this work, its derivatives are obtained directly. 

The work details the proposal and implementation of a methodology to numerically obtain first and 
second order derivatives of one variable functions. With our proposal such derivatives (without any approx- 
imation) can be obtained in straightforward manner. The methodology is easy to implement and, although 
we have coded the extended dual functions in the Fortran language, other programming languages could be 
used. 
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Table 3: Components of the acceleration vector, rgen(^) (we have taken 9 = 1 and the units are in the SI 

system) 



e 



Xgcn{9) 



fjgcnid) 



2gcn(f) 



0.00000 

0.62832 
1.25664 
1.88496 
2.51327 
3.14159 
3.76991 
4.39823 
5.02655 
5.65487 



■0.42870 
0.03898 
0.22579 
0.24420 
0.15790 
0.00803 
0.05418 
■0.02558 
0.01379 
0.14137 



0.25129 
0.15046 
0.00787 
0.11368 
0.16248 
0.19439 
0.22195 
0.14292 
0.08082 
0.34299 



0.02055 
0.56500 
0.37390 
0.23911 
0.19680 
0.09745 
0.10409 
0.31102 
0.47128 
0.56750 



Appendix A. Fortran module for the dual functions 

module sphmodual 
contains 

function rgen(th, beta, gamma, xl,x2,x3,x4) result (f_result) 
implicit none 

double precision :: f _result(3,3) 

double precision, intent(in) :: th(3) ,beta,gcimma,xl (3) ,x2(3) ,& 
x3(3) ,x4(3) 

double precision, parameter : : pi = 3 . 141592654d0 
double precision :: pids2(3) ,rcptb(3,3) ,rcptbg(3,3) 

pids2 = [pi/2d0,0d0,0d0] 

rcptb = rcp(tli,beta,xl ,x2,x3,x4) 

rcptbg = rcp(th,beta + gamma, xl ,x2,x3,x4) 

f_result = matmulvdual (rotmdual (pids2 , rcptb) , rcptbg) 

return 

end function rgen 

function rcp(th,nu,xl ,x2,x3,x4) result (f_result) 
implicit none 

double precision :: f _result(3,3) 

double precision, intent(in) :: th(3) ,nu,xl(3) ,x2(3) ,x3(3) ,x4(3) 
double precision :: nuaux(3) ,n23v(3,3) ,r2vec(3,3) 

nuaux = [nu,OdO,OdO] 

n23v = nr23(th,xl,x2,x3,x4) 

r2vec = r2(th,xl,x2) 

f_result = matmulvdual (rotmdual (nuaux, n23v) ,r2vec) 
return 

end function rep 
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function nr23(th,xl,x2,x3,x4) result (f_result) 
implicit none 

double precision :: f _result(3,3) 

double precision, intent(in) :: th(3) ,xl(3) ,x2(3) ,x3(3) ,x4(3) 
double precision :: r2v(3,3) ,r3v(3,3) ,r2cr3(3,3) ,norr2cr3(3) 

r2v = r2(th,xl,x2) 
r3v = r3(th,xl,x2,x3,x4) 
r2cr3 = crossdual(r2v,r3v) 
norr2cr3 = noniidual(r2cr3) 

f _result (1 , : ) = divisiondual(r2cr3(l , : ) ,norr2cr3) 
f_result(2, :) = divisiondual(r2cr3(2, : ) ,norr2cr3) 
f _result(3, : ) = divisiondual(r2cr3(3, : ) ,norr2cr3) 

return 

end function nr23 

function r3(th,xl ,x2,x3,x4) result (f_result) 
implicit none 

double precision :: f _result(3,3) 

double precision, intent(in) :: th(3) , xl(3), x2(3) ,x3(3) ,x4(3) 
double precision :: xld(3,3) ,x2d(3,3) ,x3d(3,3) ,x4d(3,3) ,thaux(3) 

thaux = phiodual(th,xl,x2,x3,x4) 

xld = OdO 
xld(:,l) = xl 

x2d = OdO 
x2d(: ,1) = x2 

x3d = OdO 
x3d(:,l) = x3 

x4d = OdO 
x4d(: ,1) = x4 

f_result = matmulvdual (rotmdual (thaux , x4d) , x3d) 
return 

end function r3 

function r2(th,xl,x2) result (f_result) 

implicit none 

double precision :: f _result(3,3) 

double precision, intent(in) :: th(3) , xl(3) , x2(3) 
double precision :: xld(3,3) ,x2d(3,3) 

xld = OdO 
xld(: ,l)=xl 
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x2d = OdO 
x2d(: ,1) = x2 



f_result = matmulvdual(rotmdual(th,xld) ,x2d) 
return 

end function r2 

function phiodual(th,xxl,xx2,xx3,xx4) result (f_result) 
implicit none 

double precision :: f_result(3) 

double precision, intent(in) :: th(3) ,xxl(3) ,xx2(3) ,xx3(3) ,xx4(3) 
double precision : : alphal , alpha2 , alphaS , alpha4 ,A(3),B(3),C(3), & 
auxnl2(3) ,auxn43(3) ,auxnl4(3) ,nl2(3) ,n43(3) ,nl4(3) ,& 
velmanivO) , veloscila(3) ,velfrainei(3) .velframef (3) ,angmf (3) ,angof (3) 



alphal = 


acos (Dot_Product (xxl ,xx2) ) 


alpha2 = 


acos (Dot_Product (xx2 , xx3) ) 


alph.a3 = 


acos (Dot_Product (xx3 , xx4) ) 


alpha4 = 


acos (Dot_Product (xx4 , xxl ) ) 


auxnl2 = 


cross (xxl, xx2) 


auxn43 = 


cross (xx4,xx3) 


auxnl4 = 


cross (xxl, xx4) 



nl2 = auxnl2/sqrt(suin(auxnl2**2) ) 

n43 = auxn43/sqrt (suiii(auxn43**2) ) 
nl4 = auxnl4/sqrt(suin(auxnl4**2)) 



velmaniv = velg(0d0,nl2,xxl) 
veloscila = velg(0d0,n43,xx4) 
velfrsmiei = velg(OdO,nl4,xxl) 
velfreimef = velg(alpha4,nl4,xxl) 



euigmf = [acos (Dot_Product (velmciniv, velf rgimei) ), OdO, OdO] 
angof = [acos (Dot _Product (veloscila, -velframef) ) ,OdO, OdO] 

A = sin(alphal)*sin(alph.a3)*sindual(th + angmf) 

B = [cos (alphal) *sin(alpha3) *sin(alpha4) , OdO, OdO] - sin(alphal)* & 

sin(alpha3) *cos (alpha4) *cosdual (th + angmf) 

C = sin(alphal)*cos(alpha3)*sin(alpha4)*cosdual(th + angmf) + & 
[cos (alphal) *cos (alpha3) *cos (alpha4) - cos (alpha2) , OdO , OdO] 

f_result = angof - 2dO*arctandual(divisiondual(-A-& 
sqrtdual (proddual (A , A) +proddual (B , B) -proddual (C , C) ) , C-B) ) 

return 

end function phiodual 

function cross(xx,yy) result (f_result) 
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implicit none 

double precision :: f_result(3) 

double precision, intent(in) :: xx(3) ,yy(3) 

f_result = [-(xx(3)*yy(2)) + xx(2)*yy(3) ,xx(3)*yy(l)-xx(l)*yy(3) ,& 

-(xx(2)*yy(l)) + xx(l)*yy(2)] 

return 

end function cross 

function velg (angle, vr.ra) result (f_result) 
implicit none 

double precision : : f_result(3) 

double precision, intent(in) :: angle,vr(3) ,ra(3) 

f_result=[(ra(3)*vr(2)-ra(2)*vr(3))*cos(angle)+(ra(2)*vr(l)*vr(2)+ra(3)* & 
vr(l)*vr(3) - ra(l)*(vr(2)**2+vr(3)**2))*sin(angle) , (-(ra(3)*vr(l)) + & 
ra(l)*vr(3))*cos(angle) + (vr (2) * (ra(l) *vr (1) + ra(3)*vr(3)) - ra(2)* & 
(vr(l)**2 + vr(3)**2))*sin(angle) , (ra(2)*vr(l) - ra(l)*vr(2))* & 
cos(angle)+(-(ra(3)*(vr(l)**2+vr(2)**2))+(ra(l)*vr(l)+ra(2)*vr(2))* & 
vr(3))*sin(angle) ] 

return 

end fimction velg 

function rotmdual(th, normal) result (f_result) 
implicit none 

double precision :: f _result(3,9) 

double precision, intent(in) :: th (3) , normal (3, 3) 

double precision :: nl(3) ,n2(3) ,n3(3) ,nl2(3) ,n22(3) ,n32(3) ,& 

nln2(3) ,nln3(3) ,n2n3(3) 

nl = normal ( 1 , : ) 
n2 = normal (2, :) 
n3 = normal (3 , : ) 
nl2 = proddual(nl ,nl) 
n22 = proddual(n2,n2) 
n32 = proddual(n3,n3) 
nln2 = proddual(nl,n2) 
nlnS = proddual(nl,n3) 
n2n3 = proddual(n2,n3) 

f_result(l,:) = [nl2+proddual(n22+n32,cosdual(th)) ,& 

nln2 - proddual(nln2,cosdual(th)) - proddual (n3 , sindual (th) ) ,& 
nlnS - proddual (nln3,cosdual(th)) + proddual (n2,sindual(th))] 

f _result (2 , : ) = [nln2 - proddual (nln2,cosdual(th)) + proddual (n3, s indual (th) ) ,& 
n22 + proddual (nl2+n32, cosdual (th) ) ,n2n3 -proddual(nl,s indual (th))- & 
proddual (n2n3,cosdual(th)) ] 

f _result(3, : ) = [nln3 -proddual (nln3, cosdual (th)) -proddual (n2, sindual (th) ) ,& 
n2n3+ proddual (nl , sindual (th) ) -proddual (n2n3 , cosdual (th) ) , & 
n32 + proddual (nl2+n22, cosdual (th))] 



12 



return 

end function rotmdual 

! product of a matrix with a vector (ket) 
function matmulvdual(mat,vec) result (f_re suit) 
implicit none 
!x dual 

double precision :: f _result(3,3) !y dual 
!z dual 

double precision, intent(in) :: mat(3,9), vec(3,3) 
integer : : k 



f_result = OdO 



do k=l,3 
f _result(l, 
f _result(2, 
f _result(3, 
end do 



:) = f_result(l, 
:) = f_result(2, 
:) = f_result(3, 



) + prodduaK mat(l,3*k-2:3*k) , vec(k,:) ) 
) + prodduaK mat(2,3*k-2:3*k) , vec(k,:) ) 
) + prodduaK mat(3,3*k-2:3*k) , vec(k,:) ) 



return 

end fimction matmulvdual 

function normdual(xd) result (f_result) 
implicit none 

double precision : : f_result(3) 

double precision, intent(in) :: xd(3,3) 

f_result = sqrtdual(dotdual(xd,xd)) 

return 

end function normdual 

function dotdual (xd.yd) result (f_result) 
implicit none 

double precision :: f_result(3) 

double precision, intent(in) :: xd(3,3) ,yd(3,3) 

f_result = prodduaKxdd, :) ,yd(l, :))+prodduaKxd(2, :) ,yd(2, :)) +& 
proddual(xd(3, : ) ,yd(3, : ) ) 

return 

end function dotdual 

function crossdual(xd,yd) result (f_re suit) 
implicit none 

double precision : : f _result(3,3) 

double precision, intent(in) :: xd(3,3) ,yd(3,3) 

integer : : j ,k 



f_result = OdO 
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do j=l,3 

do k=l,3 
f _result(l , : ) = f _result (1 , : ) + 
f_result(2, :) = f_result(2, : ) + 
f_result(3, :) = f_result(3, : ) + 

end do 
end do 



levic ( 1 , j , k) *proddual (xd ( j , : ) , yd (k , : ) ) 
levic (2 , j , k) *proddual (xd ( j , : ) , yd (k , : ) ) 
levic (3, j ,k)*proddual(xd(j , : ) ,yd(k, : )) 



return 

end function crossdual 

! The Levi-Civita tensor 

function levic(ii, j j ,kk) result (f_result) 

implicit none 

integer : : f _result 

integer, intent(in) :: ii,jj,kk 

integer :: j,k, auxilicir(3,3,3) ,auxiliarl(3,3) ,auxiliar2(3,3) ,& 
auxiliar3(3,3) 



auxiliarl = Reshape ( [0 ,0, , ,0, -1 ,0, 1 ,0] , [3, 3] ) 
auxiliar2 = ReshapeC [0 ,0, 1 ,0,0,0,-1 ,0,0] , [3,3] ) 
auxiliar3 = ReshapeC [0,-1 ,0, 1 ,0,0,0,0,0] , [3,3] ) 

do j=l,3 
do k=l,3 

auxiliard, j ,k) = auxiliarl (j ,k) 
end do 
end do 

do j=l,3 
do k=l,3 

auxiliar(2, j ,k) = auxiliar2(j ,k) 
end do 
end do 

do j=l,3 
do k=l,3 

auxiliar(3, j ,k) = auxiliar3(j ,k) 
end do 
end do 

f_result = auxiliax(ii, j j ,kk) 
return 

end function levic 

function divisiondual(x,y) result (f_result) 
implicit none 

double precision :: f_result(3) 

double precision, intent(in) :: x(3),y(3) 
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f_result=[ x(l)/y(l),( y(l)*x(2)-x(l)*y(2) )/(y(l)**2) ,& 

x(3)/y(l)-2d0*x(2)*y(2)/y(l)**2d0+x(l)*(2d0*y(2)**2/y(l)**3d0- y(3)/y(l)**2d0) ] 



return 

end function divisiondual 

function sqrtdual(x) result (f_result) 
implicit none 

double precision : : f_result(3) 

double precision, parameter :: mediodr(3)=[0.5d0,0d0,0d0] 
double precision, intent(in) :: x(3) 
f _result=powdual (x.mediodr) 
return 

end function sqrtdual 

function powdual(x,y) result (f_result) 
implicit none 

double precision :: f _result (3) ,ml ,m2 
double precision, intent(in) :: x(3),y(3) 
double complex :: xz, yz, zl, z2, z3 

double precision, parameter : : eps = ld-12 

! depending on the machine precision, the eps number can be made smaller 

xz = x(l) + (OdO,OdO) 
yz = yd) + (OdO,OdO) 

7.1= x(l)**y(l) 

z2 = zl*(x(2)*y(l)/x(l) + y(2)*log( xz ) ) 

z3 = zl*((x(2)**2*(-ld0 + y(l))*y(l))/x(l)**2 + (x(3)*y(l)+ & 
2d0*x(2)*(ld0 + Log(xz)*y(l))*y(2))/x(l) + Log(xz)*(Log(xz)* & 
y(2)**2 + y(3))) 

ml=OdO 
m2=0d0 

f_result = ml/m2 

! dropping small complex parts 

if ( abs (aimag(zl) ) . le . eps) f_result(l) = dble(zl) 
if( abs(aimag(z2)) .le.eps) f_result(2) = dble(z2) 
if( abs (aimag(z3)) .le.eps) f_result(3) = dble(z3) 

return 

end function powdual 

function proddual(x,y) result (f_result) 
implicit none 

double precision :: f_result(3) 

double precision, intent(in) :: x(3) , y(3) 

f_result=[x(l)*y(l) , x(2)*y(l) + x(l)*y(2), x(3)*y(l)+2d0*x(2)* & 

y(2)+x(l)*y(3)] 

return 

end function proddual 
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function sindual(x) result (f_re suit) 
implicit none 

double precision :: f_result(3) 
double precision, intent(in) :: x(3) 

f_result=[sin(x(l)) ,cos(x(l))*x(2) ,-sin(x(l))*x(2)**2 + & 

cos(x(l))*x(3)] 

return 

end function sindual 

function cosdual(x) result (f_re suit) 
implicit none 

double precision :: f_result(3) 
double precision, intent(in) :: x(3) 

f_result=[cos(x(l)) ,-sin(x(l))*x(2) ,-cos(x(l))*x(2)**2 -& 

sin(x(l))*x(3)] 

return 

end function cosdual 

function tcindual(x) result (f_re suit) 

implicit none 

double precision :: f_result(3) 
double precision, intent(in) :: x(3) 

f_result=[tan(x(l)) ,x(2)/cos(x(l))**2, x(3)/cos(x(l))**2 + & 

(2d0*tan(x(l))/cos(x(l))**2)*x(2)**2] 

return 

end function tcmdual 

function arcsindual(x) result (f_result) 

implicit none 

double precision :: f_result(3) 
double precision, intent(in) :: x(3) 

f_result=[asin(x(l)) ,x(2)/sqrt(l-x(l)**2) ,x(3)/sqrt(l-x(l)**2) + & 

(x(l)/(l-x(l)**2)**(1.5d0))*x(2)**2] 

return 

end function arcsindual 

function arccosdual(x) result (f_result) 

implicit none 

double precision :: f_result(3) 
double precision, intent(in) :: x(3) 

f _result=[acos(x(l)) ,-x(2)/sqrt (l-x(l)**2) ,-x(3)/sqrt(l-x(l)**2)-& 

(x(l)/(l-x(l)**2)**(1.5d0))*x(2)**2] 
return 

end function arccosdual 

function arctandual(x) result (f_result) 

implicit none 

double precision :: f_result(3) 
double precision, intent(in) :: x(3) 

f_result=[atan(x(l)) ,x(2)/(ld0+x(l)**2) , x(3)/(ld0+x(l)**2) + & 
(-2d0*x(l)/(x(l)**2+l)**2)*x(2)**2 ] 
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return 

end function arctandual 

function expdual(x) result (f_result) 
implicit none 

double precision :: f_result(3) 
double precision, intent(in) :: x(3) 

f_result=[exp(x(l)) ,exp(x(l))*x(2) ,exp(x(l))*x(2)**2 + exp(x(l))*& 

x(3) ] 

return 

end function expdual 

function logdual(x) result (f_result) 
implicit none 

double precision :: f_result(3) 
double precision, intent(in) :: x(3) 

f_result=[log(x(l)) , x(2)/x(l), x(3)/x(l) - x(2)**2/x(l)**2] 
return 

end function logdual 

function erf dual (x) result (f_re suit) 
implicit none 

double precision :: f_result(3) 

double precision, intent(in) :: x(3) 

f_result=[erf (x(l)) ,2d0/sqrt(2d0*asin(ld0))*exp(-x(l)**2)*x(2) ,& 
x(3) *2d0/sqrt (2d0*asin(ld0) ) *exp(-x(l) **2) -& 
4dO*x (1) /sqrt (2d0*asin(ld0) ) *exp(-x(l) **2) *x(2) **2] 
return 

end function erfdual 

function shdual(x) result (f_result) 
implicit none 

double precision :: f_result(3) 
double precision, intent(in) :: x(3) 

f_result=[sinh(x(l)) , cosh(x(l) ) *x(2) , x(3)*cosh(x(l)) +& 

sinh(x(l))*x(2)**2] 

return 

end fimction shdual 

function chdual(x) result (f_result) 
implicit none 

double precision :: f_result(3) 
double precision, intent(in) :: x(3) 

f_result=[cosh(x(l)) , sinh(x(l) ) *x(2) , x(3)*sinli(x(l)) +& 

cosh(x(l))*x(2)**2] 

return 

end function chdual 

function tanhdual(x) result (f_result) 
implicit none 

double precision :: f_result(3) 
double precision, intent(in) :: x(3) 
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f_result=[tanh(x(l)) , x(2)/cosh(x(l))**2,x(3)/cosh(x(l))**2 + & 

(-2d0*x(2)**2)*tanh(x(l))/cosh(x(l))**2 ] 

return 

end fimction tanhdual 

end module sphmodual 



Appendix B. User guide 

Here we briefly describe the functions of the module. The kind of variables can be seen in the code of 
the module. For the functions related to the spherical mechanism, a unit sphere is assumed. 

• tanhdual (x) 

dual version for the Tangent function. 

• chdual(x) 

dual version for the Hyperbolic cosine function. 

• shdual(x) 

dual version for the Hyperbolic sine function. 

• erf dual (x) 

dual version for the error function. 

• logdual(x) 

dual version for the logarithmic function. 

• expdual(x) 

dual version for the exponential function. 

• arctandual(x) 

dual version for the arctangent function. 

• arccosdual (x) 

dual version for the arccosine function. 

• arcsindual (x) 

dual version for the arcsine function. 

• tandual(x) 

dual version for the tangent function. 

• cosdual(x) 

dual version for the cosine function. 

• sindual(x) 

dual version for the sine function. 

• proddual(x,y) 

dual version for the product of two dual functions. 

• powdual(x,y) 

dual version for the power of two dual functions (a;^). 

• sqrtdual(x) 

dual version for the square root function. 
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• divisiondual(x,y) 

dual version for the division of two dual functions (x/y). 

• matmulvdual(mat,vec) 

dual version for the multiplication of a matrix with a vector. The matmulvdual function is an array 
of dimension 3x3. 

Notice that this is not in general the matrix product, but the product of a matrix with a vector (by 
vector we mean a 3D-vector not the abstract definition of a vector) . 

• crossdual (xd,yd) 

dual version for the cross product. 

• dotdual(xd,yd) 

dual version for the dot product. 

• normdual(xd) 

dual version for the euclidean norm of a vector. 

• rotmdual(th,uaxis) 

dual version for the rotation matrix of dual angle th in direction of the dual unit vector uaxis. 

• levied, j ,k) 

The Levi-Civitta tensor. This is not a dual function but is required for some of them. 

• velg (angle, vr,ra) 

The first derivative of the parametric equation for the geodesic on the unit sphere connecting the 
points vr and va. This is not a dual function. It is useful to calculate angles between geodesies. 

• cross(x,y) 

Cross product of the 3D-vectors x and y. This is not a dual function, is the common cross product. 

• phiodual(thd,xdl,xd2,xd3,xd4) 

dual version of the output angle (Eq.(15) of the paper). 

• r2(th,xl,x2) 

dual version for the unit position vector of the input link (Eq. (12) of the paper). 

• r3(th,xl,x2,x3,x4) 

dual version for the unit position vector of the output link (Eq. (13) of the paper). 

• nr23(th,xl,x2,x3,x4) 

dual version for the unit vector normal to r2 and (Eq. (20) of the paper). 

• rcp(th,nu,xl,x2,x3,x4) 

dual version of the Eq. (19) of the paper. 

• rgen(th, beta, gamma, xl,x2,x3,x4) 

dual version of the generated coupler-point curve. Eq. (21) of the paper. 

Below a simple Fortran program to calculate the first and second derivative dX x — 0.5 of the cross 
product of the vectors v(a::) — {erfix) sin(a::), cos(a;) + x"^,e^} and w{x) — {x — cos(a;), sin(2;), 1/x} is 
presented. The error function is defined as: 
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program crossp 
use sphmodual 
implicit none 

double precision :: ang(3) ,xcl(3,3) , xc2(3,3) ,thr(3) ,one(3) ,crsp(3,3) 

thr = [3d0,0d0,0d0] ! dual version for the number 3 
one = [Id0,0d0,0d0] ! dual version for the number 1 
ang = [0.5d0,ld0,0d0] 



xcKl, 


) 


= proddual( erf dual (ang) ,sindual(axLg)) 


xcl(2. 


) 


= cosdualCcoig) + powdualCcOig.thr) 


xcl(3. 


) 


= expdual(ang) 


xc2(l. 


) 


= ang - cosdual(ang) 


xc2(2, 


) 


= sindual(cing) 


xc2(3. 


) 


= divisiondual(one,ang) 



crsp = crossdual(xcl,xc2) 



print*, "The cross product evaluated in 0.5 is:" 

print*, crsp(:,l) 

print*, "The first derivative evaluated in 0.5 is:" 

print*, crsp(:,2) 

print*, "The second derivative evaluated in 0.5 is:" 

print*, crsp(:,3) 

end program crossp 

Saving the module and the above program in the files sphmodual. f 90 and crossp. f 90 respectively, we 
can compile, using the gfortran compiler, as follows: gfortran sphmodual. f 90 crossp. f 90 
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